Day07:地球観測・科学データの可視化
はじめに
本日は地震観測データとその他の観測データを組合せて可視化する手法を学びます．
地震観測データには，アメリカ地質調査所(United States Geological Survey; USGS) Earthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)が公開するデータを用います．
USGSの地震観測データと組み合わせるその他の観測データには， テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データとスミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM(http://www.volcano.si.edu/)が公開する火山データを用います．
本日の演習では，緯度（latitude）・経度（longitude）という情報を基に複数のデータを組み合わせる手法の習得を目標とします．
プレート境界線の可視化
まず，テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データを可視化します。
プレート境界線とはプレートテクトニクス（プレート理論）により説明されるプレートとプレートの境目のことです．詳しくは以下リソースを参照してください．
http://ja.wikipedia.org/wiki/プレートテクトニクス
プレート境界線データのダウンロード
テキサス大学地球物理研究所PLATESプロジェクトのデータ公開ページ（http://www-udc.ig.utexas.edu/external/plates/data.htm）からプレート境界線データをダウンロードします．
ページの中頃に「PLATES_PlateBoundary_ArcGIS.zip」というファイルがあるので，それをダウンロードします。
Rの作業ディレクトリにディレクトリ「plate」を作成し， そのディレクトリに「PLATES_PlateBoundary_ArcGIS.zip」を展開してください．
Rで上記処理をさせる場合は以下となります．
# fileがなければダウンロードする
if(!file.exists("PLATES_PlateBoundary_ArcGIS.zip")){
download.file(url="http://www-udc.ig.utexas.edu/external/plates/data/plate_boundaries/PLATES_PlateBoundary_ArcGIS.zip",destfile = "PLATES_PlateBoundary_ArcGIS.zip")
}
# directoryがなければファイルを展開する
if(!dir.exists("plate")){
unzip("PLATES_PlateBoundary_ArcGIS.zip",exdir = "plate")
}ダウンロードしたプレート境界線データはESRI Shapeフォーマットで格納されています． Shape形式はGIS(地理情報システム)で比較的よく利用されているファイルフォーマットです． 詳しくは Esri Japanのシェープファイルについて等を参照してください．
Shapeファイル処理用ライブラリのインストール
Shapeファイルを処理するため maptools ライブラリをインストールします． また，データフレーム変換用のライブラリ broomもインストールします．
#ライブラリの読み込み．なければインストールする．
if(!require(ggplot2)){
install.packages("ggplot2")
library(ggplot2)
}
if(!require(maps)){
install.packages("maps")
library(maps)
}
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
if(!require(broom)){
install.packages("broom")
library(broom)
}
#游ゴシック体を使う
if(.Platform$OS.type=="windows")
windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))Shapeファイルの描画
それでは早速Shapeファイルを描画しましょう．
発散型境界の描画
展開したShapeファイルを読み込み，ggplotで発散型境界を描画します．
ggplot2 は data.frame でないと処理できないので， maptoolsの処理結果（SpatialLinesDataFrame）をそのまま利用することができません． そこで，broom パッケージのtidy 関数を利用して SpatialLinesDataFrame を data.frame に変換ます．
それでは早速下記コードを実行してください．
# shapeファイルの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
# オブジェクトのクラスを確認する
class(ridge)## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"ridge.fort<-broom::tidy(ridge)
class(ridge.fort)## [1] "data.frame"g <- ggplot() + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("発散型境界 (data by UTIG Plate Proj.)")
g
発散型境界線が描画されました．
世界地図と合わせて描画してみましょう．
g+borders("world")
トランスフォーム型境界の描画
つぎにトランスフォーム型境界を描画します。
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
g <- ggplot() + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("トランスフォーム型境界 (data by UTIG Plate Proj.)")
g
これも世界地図と合わせて描画してみましょう．
g+borders("world")
収束型境界の描画
つぎに収束型境界を描画します。
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
g <- ggplot() + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("収束型境界 (data by UTIG Plate Proj.)")
g
収束型境界についても，世界地図と合わせて描画してみましょう．
g+borders("world")
すべての境界線の描画
それでは発散型，トランスフォーム型，収束型境界線を合わせて描画します．
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
g <- ggplot()
g <- g + theme_bw(base_family = "yugo")+ggtitle("プレート境界 (data by UTIG Plate Proj.)")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
g
世界地図と合わせます．
g + borders(database = "world")
地震はプレート境界線で生じる歪が原因となり発生するといわれています．
日本で地震が多い理由も納得できますね．
地震観測データとプレート境界線の組合せ
地震観測データの描画
アメリカ地質調査所(United States Geological Survey; USGS) のEarthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)の地震観測データとテキサス大学地球物理研究所PLATESプロジェクトのプレート境界データを合成します．
まずは地震観測データを描画します． データは前回の可視化で用いた「20110311JSTquery.csv」を使います．
csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
#地震をプロット（色で深度，サイズでマグニチュードを表現）
g<-ggplot(csv,aes(x=longitude,y=latitude,color=depth,size=mag))
g<-g+geom_point()
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ（data by USGS Earthquake Hazards Prog. Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く，深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
g
ここまでは前回と同じですね．
地震観測データとプレート境界線データの描画
それでは地震観測データとプレート境界線データと合わせて可視化してみましょう．
ポイントはこれまで描画データを指定していたggplot()関数では何も指定せずに， 点や線を描画する geom_point や geom_path で描画データを指定している点です． こうすることで複数のデータを用いた可視化が可能となります（ggplot 関数でデータを指定すると全体でデータが共有されることになる）．
# 地震観測データの読み込み
csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
# プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
# 描画データを指定しない
g<-ggplot()
# geom_point で描画するデータ（地震観測データ）を指定している．
g<-g+geom_point(data = csv,mapping = aes(x=longitude,y=latitude,color=depth,size=mag))
# 世界地図の描画
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ（data by USGS Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く，深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
# タイトルを上書き
g <- g +ggtitle("プレート境界と地震データ(data by UTIG Plate Proj. and USGS)")
#発散型境界の描画（データを指定している点に注目）
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画（データを指定している点に注目）
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画（データを指定している点に注目）
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
g
プレート境界線近辺で地震が発生していることがわかります．
地震発生メカニズム（気象庁，地震発生のしくみ）を実際のデータで確認することができました．
プレート境界線データと火山データの組合せ
つぎにプレート境界線データと火山データを組合せて可視化してみます．
火山もプレートテクトニクスに関連が強い地形です．（参考：Wikipdia，火山）
火山データには，スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM（http://www.volcano.si.edu/）が公開する火山データを用います．
火山データのダウンロードとフォーマット変換
火山データのダウンロードするため，スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMサイト(http://www.volcano.si.edu/)へアクセスします．
「Database」 > 「HOLOCENE SPREADSHEET」をクリックします．
「GVP_Volcano_List.xls」というファイル名でダウンロードが始まります．
ダウンロードが完了したら，「GVP_Volcano_List.xls」をCSVフォーマットに変換します．
フォーマット変換にはExcelを利用します．Excelでファイルを開き，「名前を付けて保存」でCSVフォーマットとして保存します．
変換後のCSVファイル「GVP_Volcano_List.csv」をRの作業ディレクトリに移動します．
火山データの描画
それではさっそく火山データを読み込んでみましょう．
#データを読み込む
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1,stringsAsFactors = FALSE)
#構造を確認する
str(volc)## 'data.frame': 1520 obs. of 13 variables:
## $ Volcano.Number : int 210010 210020 210030 210040 211001 211003 211004 211010 211020 211030 ...
## $ Volcano.Name : chr "West Eifel Volcanic Field" "Chaine des Puys" "Olot Volcanic Field" "Calatrava Volcanic Field" ...
## $ Country : chr "Germany" "France" "Spain" "Spain" ...
## $ Primary.Volcano.Type: chr "Maar(s)" "Lava dome(s)" "Pyroclastic cone(s)" "Pyroclastic cone(s)" ...
## $ Activity.Evidence : chr "Eruption Dated" "Eruption Dated" "Evidence Credible" "Eruption Dated" ...
## $ Last.Known.Eruption : chr "8300 BCE" "4040 BCE" "Unknown" "3600 BCE" ...
## $ Region : chr "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" ...
## $ Subregion : chr "Western Europe" "Western Europe" "Western Europe" "Western Europe" ...
## $ Latitude : num 50.2 45.8 42.2 38.9 43.2 ...
## $ Longitude : num 6.85 2.97 2.53 -4.02 10.87 ...
## $ Elevation..m. : int 600 1464 893 1117 500 800 949 458 1281 789 ...
## $ Dominant.Rock.Type : chr "Foidite" "Basalt / Picro-Basalt" "Trachybasalt / Tephrite Basanite" "Basalt / Picro-Basalt" ...
## $ Tectonic.Setting : chr "Rift zone / Continental crust (>25 km)" "Rift zone / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" ...#要約を確認する
summary(volc)## Volcano.Number Volcano.Name Country
## Min. :210010 Length:1520 Length:1520
## 1st Qu.:261148 Class :character Class :character
## Median :300020 Mode :character Mode :character
## Mean :296768
## 3rd Qu.:342123
## Max. :600000
## Primary.Volcano.Type Activity.Evidence Last.Known.Eruption
## Length:1520 Length:1520 Length:1520
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## Region Subregion Latitude Longitude
## Length:1520 Length:1520 Min. :-78.50 Min. :-179.97
## Class :character Class :character 1st Qu.: -7.13 1st Qu.: -78.15
## Mode :character Mode :character Median : 13.84 Median : 38.39
## Mean : 14.05 Mean : 23.85
## 3rd Qu.: 42.00 3rd Qu.: 139.04
## Max. : 85.61 Max. : 179.58
## Elevation..m. Dominant.Rock.Type Tectonic.Setting
## Min. :-4200.0 Length:1520 Length:1520
## 1st Qu.: 680.8 Class :character Class :character
## Median : 1462.0 Mode :character Mode :character
## Mean : 1677.0
## 3rd Qu.: 2344.8
## Max. : 6879.0#データを確認する
head(volc)## Volcano.Number Volcano.Name Country Primary.Volcano.Type
## 1 210010 West Eifel Volcanic Field Germany Maar(s)
## 2 210020 Chaine des Puys France Lava dome(s)
## 3 210030 Olot Volcanic Field Spain Pyroclastic cone(s)
## 4 210040 Calatrava Volcanic Field Spain Pyroclastic cone(s)
## 5 211001 Larderello Italy Explosion crater(s)
## 6 211003 Vulsini Italy Caldera
## Activity.Evidence Last.Known.Eruption Region
## 1 Eruption Dated 8300 BCE Mediterranean and Western Asia
## 2 Eruption Dated 4040 BCE Mediterranean and Western Asia
## 3 Evidence Credible Unknown Mediterranean and Western Asia
## 4 Eruption Dated 3600 BCE Mediterranean and Western Asia
## 5 Eruption Observed 1282 CE Mediterranean and Western Asia
## 6 Eruption Observed 104 BCE Mediterranean and Western Asia
## Subregion Latitude Longitude Elevation..m.
## 1 Western Europe 50.170 6.85 600
## 2 Western Europe 45.775 2.97 1464
## 3 Western Europe 42.170 2.53 893
## 4 Western Europe 38.870 -4.02 1117
## 5 Italy 43.250 10.87 500
## 6 Italy 42.600 11.93 800
## Dominant.Rock.Type
## 1 Foidite
## 2 Basalt / Picro-Basalt
## 3 Trachybasalt / Tephrite Basanite
## 4 Basalt / Picro-Basalt
## 5 No Data (checked)
## 6 Trachyte / Trachydacite
## Tectonic.Setting
## 1 Rift zone / Continental crust (>25 km)
## 2 Rift zone / Continental crust (>25 km)
## 3 Intraplate / Continental crust (>25 km)
## 4 Intraplate / Continental crust (>25 km)
## 5 Subduction zone / Continental crust (>25 km)
## 6 Subduction zone / Continental crust (>25 km)#富士山のデータを確認する
volc[volc$Volcano.Name=="Fujisan",]## Volcano.Number Volcano.Name Country Primary.Volcano.Type
## 598 283030 Fujisan Japan Stratovolcano
## Activity.Evidence Last.Known.Eruption Region
## 598 Eruption Observed 1708 CE Japan, Taiwan, Marianas
## Subregion Latitude Longitude Elevation..m. Dominant.Rock.Type
## 598 Honshu 35.361 138.728 3776 Basalt / Picro-Basalt
## Tectonic.Setting
## 598 Subduction zone / Continental crust (>25 km)火山データの位置が緯度（Latitude)，経度（Longitude）に記録されていることがわかります．
それでは火山データを可視化してみましょう．ここでは散布図 geom_point を用いていますが，単なる点だと火山ぽくないので，「shape=24」を指定することで点を三角に変更しています．
#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
geom_point で指定可能な shape には他にも以下のような形があります．
g<-ggplot(data=data.frame(x=0:24))+scale_shape_identity()
g<-g+geom_point(mapping=aes(x=x,y=0,shape=x),size=5)
g<-g+geom_point(mapping=aes(x=x,y=1,shape=x),size=5,fill="red")
g<-g+xlab("shape types")
g
今回はshapeの24番（色塗り可能な三角形）を指定することで火山を表現しています．
スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMの火山データにはその他にも様々な情報が記録されています．例えば，「Elevation..m.」には火山の標高が記録されています．
可視化してみると次のようになります．
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude,color=Elevation..m.,fill=Elevation..m.),shape=24)+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program)")
火山データとプレート境界線データの描画
それでは火山データとプレート境界線データを合わせて描画してみましょう．
まず，火山データと世界地図を描画します．
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
g<-ggplot()
g<-g+borders(database = "world")
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
g
この描画結果にプレート境界線を追加すればよいわけです．
地震観測データとプレート境界線データの描画を参考にすると， 下記のようにコードを書くことができると思います．要は描画単位（geom_pointやgeom_pathごと）に可視化するデータを指定すればよいのです．
#プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
#火山データの読み込み
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
g<-ggplot()
#地図の描画
g<-g+borders(database = "world")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
#火山の描画
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
g
このように火山もプレート境界線近辺に多くできていることがわかります．
まとめ
本日は地球観測データと科学データを複数組合せて可視化する手法を学びました． 具体的には，緯度（latitude）・経度（longitude）をキーに，地震観測データ，プレート境界線データ，火山データを関連付けて可視化しました．このように複数のデータを組合せて可視化するにはデータを関連付けるためのキーが必要となります．今回の演習で感じてもらえたと思いますが，緯度・経度は複数データを組み合わせる際のキー候補としやすい情報の一つです．
次回は今回の内容に加えて地理空間情報の可視化に取り組みます．
課題
以下の問題を解いてください． その際に用いたRスクリプトファイル，結果として得られた画像ファイルをZIPでまとめて提出してください．
※Rスクリプトファイルは可視化毎に作成してください．可視化に必要な内容のみを含めるようにしてください． ※データファイルは不要です．
【問題】地震，プレート境界線，火山データのすべてを可視化してください
ヒント：
- 火山データの標高は可視化しなくて構いません
- 見やすくなるよう描画順や色等を工夫してください
可視化例：
【発展問題】上記内容を特定の国に絞って可視化してください
ヒント：
- ggplot2 の可視化データはすべてdata.frameです
- subsetを使います
例：
Masaharu Hayashi を著作者とするこの 作品 は クリエイティブ・コモンズの 表示 4.0 国際 ライセンスで提供されています。